clc;clear;close;
% 添加必要库
addpath(genpath('../utils'));

%% 参数设置
% 原图
cb_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\cb.png"; % 棋盘格原图
noise_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\noise1.png"; % 噪声原图
black_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\black.png"; % 黑色原图
white_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\white.png"; % 白色原图
% 拍摄图
cb = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\cb\光谱图像\550nm(3.65mm).tif";% 拍摄的棋盘格
noise = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\noise1\550nm聚集位置所有中心波长图像\650nm(3.65mm).tif"; % 拍摄的噪声图像
black = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\black\550nm聚集位置所有中心波长图像\650nm(3.65mm).tif"; % 拍摄的黑板
white = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\white\550nm聚集位置所有中心波长图像\650nm(3.65mm).tif";% 拍摄的白板
% 
R = 21; % 模糊半径
cb_rows = 15;
cb_cols = 20;
rect = [288,197,154,105];

%% 1.读入原图与拍摄图
% 原图
cb_prime = im2double(imread(cb_prime));

noise_prime = im2double(imread(noise_prime)); % 拍摄的噪声图像
noise_prime = imresize(noise_prime,size(cb_prime),'nearest');% 将噪声图像缩放到与cb_prime同样大小

black_prime = im2double(imread(black_prime)); % 拍摄的黑板
black_prime = imresize(black_prime,size(cb_prime),'nearest');% 将噪声图像缩放到与cb_prime同样大小

white_prime = im2double(imread(white_prime));% 拍摄的白板
white_prime = imresize(white_prime,size(cb_prime),'nearest');% 将噪声图像缩放到与cb_prime同样大小

% 聚焦位置数据
img_cb = im2double(imread(cb));

% 离焦位置数据
img_noise = im2double(imread(noise));
img_blk = im2double(imread(black));
img_wht = im2double(imread(white));

figure('units','normalized','outerposition',[0 0 1 1])
subplot(4,2,1)
imshow(cb_prime)
title('棋盘格原图');
subplot(4,2,3)
imshow(noise_prime)
title('伯努利噪声原图');
subplot(4,2,5)
imshow(black_prime)
title('黑板原图');
subplot(4,2,7)
imshow(white_prime)
title('白板原图');
subplot(4,2,2)
imshow(img_cb)
title('棋盘格拍摄图');
subplot(4,2,4)
imshow(img_noise)
title('伯努利噪声拍摄图');
subplot(4,2,6)
imshow(img_blk)
title('黑板拍摄图');
subplot(4,2,8)
imshow(img_wht)
title('白板原图');
pause(2);
close;

%% 3.进行图像配准
% iprime_norm = warp_function(cb_prime,...
%                 noise_prime,...
%                 img_cb,...
%                 cb_rows,...
%                 cb_cols);
iprime_norm = img_cb; % 直接使用

%% 4.进行明暗校正
u_img = img_blk + iprime_norm .* max((img_wht - img_blk),0);

figure('units','normalized','outerposition',[0 0 1 1])
subplot(1,2,1)
imshow(u_img )
title('潜在清晰图像');
subplot(1,2,2)
imshow(img_noise)
title('拍摄图像');
pause(2)
close;
% save('results/subset_20220113_uimgs_bimgs.mat','img_noise','u_img')

%% 5.进行PSF估计
k_size = 3*R;

% Load the u_imgs and b_imgs for PSF estimation
% load(fullfile('results','subset_20220113_uimgs_bimgs.mat'));
b_img = img_noise;
u_img = u_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);
b_img = b_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);

% k_init = fspecial('gaussian',k_size,0.5);
k_init = zeros(k_size); k_init(ceil(k_size/2), ceil(k_size/2)) = 1;   
tic
psf = PSF_estimation(k_init,u_img,b_img,R);
toc

%% 6. 对PSF进行边缘零填充
psf_size = 101;
psf = padarray(psf,[(psf_size-k_size)/2,(psf_size-k_size)/2],'both');
figure('units','normalized','outerposition',[0 0 1 1])
imshow(psf,[])
exportgraphics(gcf,'PSF.png')


%% 7.进行PSF去模糊测试结果
img = imread("raw_data\20220113小测试\fen 550.tif");
img = double(img)/4095;
img = imcrop(img,rectangleROI);
img_deblur = deconvlucy(img,psf);
figure()
subplot(1,2,1)
imshow(img);
title('拍摄图');
subplot(1,2,2)
imshow(img_deblur);
title('去模糊图');

% 后续可删的图像输出函数
imwrite(u_img,'u_img.png');
imwrite(b_img,'b_img.png');